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Monte Carlo simulations and finite-size scaling analysis have been carried out to study the critical 
behavior in a two-dimensional system of particles with two bonding sites that, by decreasing temper- 
ature or increasing density, polymerize reversibly into chains with discrete orientational degrees of 
freedom and, at the same time, undergo a continuous isotropic-nematic (IN) transition. A complete 
phase diagram was obtained as a function of temperature and density. The numerical results were 
compared with Mean Field (MF) and Real Space Renormalization Group (RSRG) analytical pre- 
dictions about the IN transformation. While the RSRG approach supports the continuous nature 
of the transition, the MF solution predicts a first-order transition line and a tricritical point, at 
variance with the simulation results. 

PACS numbers: 05.50.+q, 64.70. mf, 61.20.Ja, 64.75.Yz, 75.40.Mg 



I. INTRODUCTION 

Molecular self-assembly is one of the basic mechanisms 
of life and matter, and thus, modeling and measurements 
of naturally occurring self-assembling systems has long 
been pursued in the biological and physical sciences [l|, ■ 
Despite the large number of papers that are currently re- 
ported, many of the ideas that are crucial to the devel- 
opment of this area (molecular shape, interplay between 
enthalpy and entropy, nature of the forces that connect 
the particles in self-assembled molecular aggregates) are 
simply not yet under the control of investigators. 

Self-assembly also poses a number of substantial tech- 
nological challenges |5T-fl3|. In fact, the biological sys- 
tems use self-assembly to assemble macromolecules and 
structures. Imitating these strategies and creating 
novel molecules with the ability to self-assemble into 
supramolecular assemblies is an important technique in 
nanotcchnology. There is then a need for understanding 
the basic principles governing this type of organization. 

It is obvious that a complete analysis of the self- 
assembly phenomenon is a quite difficult subject because 
of the complexity of the involved microscopic mecha- 
nisms. For this reason, the understanding of simple mod- 
els with increasing complexity might be a help and a 
guide to establish a general framework for the study of 
this kind of systems, and to stimulate the development 
of more sophisticated models which can be able to repro- 
duce concrete experimental situations. 

Computer simulations have shown that spherical par- 
ticles interacting isotropically through repulsive inter- 
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particle interactions can spontaneously assemble into 
anisotropic structures The presence of an 

isotropic short ranged interparticle attraction coupled to 
a longer ranged repulsion can also yield anisotropic struc- 
tures. However, most real components, from proteins 
to ions [ljj to the wide variety of recently synthesized 
nanoparticles [H, [s| , interact via anisotropic or "patchy" 
attractions. Simulation work [l8l - l22j reveals assembly 
pathways of such components to be in general richer than 
those of their isotropic counterparts. Experimental real- 
ization of such systems is growing. An example of real 
patchy particles is presented in Ref. 23] . Such particles 
offer the possibility to be used as building blocks of specif- 
ically designed self-assembled structures [1, H, d, Hj, [25[ . 
Moreover, the implications of patchy colloids for proteins 
plj . which are patchy by nature, could be significant. 

In this line, we consider in this paper the general 
problem of particles with strongly anisotropic, highly 
directional interactions in which effectively attractive 
patches induce the reversible self-assembly of particles 
into chains, i.e., equilibrium polymerization [26l432| . Re- 
cently, several research groups reported on the assembly 
of colloidal particles in linear chains. Selectively function- 
alizing the ends of hydrophilic nanorods with hydropho- 
bic polymers, Nie et al. reported the observation of rings, 
bundles, chains, and bundled chains (33[. In another 
experimental study carried out by Chang et al. [34j ]. 
gold nanorods were assembled into linear chains using 
a biomolecular recognition system. In a direct relation 
with the present work, Clair et al. [35| investigated the 
self-assembly of terephthalic acid (TPA) molecules on the 
Au(lll) surface. Using scanning tunneling microscopy, 
the authors showed that the TPA molecules arrange in 
one-dimensional chains with a discrete number of orien- 
tations relative to the substrate. 
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FIG. 1: Schematic representation of a system of self- 
assembled rigid rods on a square lattice. 

It is well known that solutions of self-assembled chains 
exhibit a transition from a disordered isotropic phase to 
an ordered nematic phase as the concentration of par- 
ticles increases. Experimental examples of equilibrium 
polymer systems that exhibit a isotropic-ncmatic (IN) 
phase transition include wormlike micelles 13611 and self- 
assembled protein fibers like /-actin [§3, l38j . In this 
context, a recent paper was devoted to the study of a 
system of self-assembled rigid rods adsorbed on a two- 
dimensional lattice [13] ■ In Ref. [13], Tavares et al. 
studied a system composed of monomers with two attrac- 
tive (sticky) poles that polymerize reversibly into polydis- 
perse chains and, at the same time, undergo a isotropic- 
nematic (IN) continuous phase transition [39j]. So, the 
interplay between the self-assembly process and the ne- 
matic ordering is a distinctive characteristic of these sys- 
tems. Using an approach in the spirit of the Zwanzig 
model (45j , the authors found that nematic ordering en- 
hances bonding. In addition, the average rod length was 
described quantitatively in both phases, while the loca- 
tion of the ordering transition, which was found to be 
continuous, was predicted semiquantitatively by the the- 
ory. With respect to the characteristics of the phase tran- 
sition, it has recently been shown that, at intermediate 
density, the IN transition is in the q = 1 Potts universal- 
ity class [11]. 

The temperature-coverage phase diagram obtained in 
Ref. [13] is qualitative only, and the theory overestimates 
the critical temperature in all range of coverage. In addi- 
tion, the possibility of a reentrant nematic transition at 
high densities [46| was not investigated by Tavares et al.. 
Accordingly, the main objective of the present work is to 



provide an accurate determination of the phase diagram 
of the system. For this purpose, extensive Monte Carlo 
(MC) simulations, supplemented by finite-size scaling 
analysis and two analytical approximations, have been 
carried out to obtain the critical temperature character- 
izing the IN phase as a function of the coverage. The 
paper is organized as follows. In Sec. II we describe the 
lattice-gas model. The simulation scheme and computa- 
tional results are given in Sec. III. In Sec. IV we present 
the analytical approximations [Mean-Field approxima- 
tion (MF), and Real Space Renormalization Group ap- 
proach (RSRG)] and compare the MC results with the 
theoretical calculations. Finally, the general conclusions 
are drawn in Sec. V. 



II. LATTICE-GAS MODEL 

As in Refs. [13, Hil, we consider a system of self- 
assembled rods with a discrete number of orientations 
in two dimensions. We assume that the substrate is rep- 
resented by a square lattice of M = L x L adsorption 
sites, with periodic boundary conditions. N particles 
are adsorbed on the substrate with two possible orienta- 
tions along the principal axis of the square lattice. These 
particles interact with nearest-neighbors (NN) through 
anisotropic attractive interactions (see Fig. 1). Then, 
the adsorbed phase is characterized by the Hamiltonian 

H = -w^2 Vii ■ °i ■ (1) 

where indicates a sum over NN sites; w represents 
the NN lateral interaction between two neighboring i and 
j, which are aligned with each other and with the inter- 
molecular vector fij ; and &i is the occupation vector with 
&i = if the site i is empty, c?; = x if the site i is oc- 
cupied by a particle with orientation along the x-axis, 
and c?j = y if the site i is occupied by a particle with 
orientation along the y-axSs. 

A cluster or uninterrupted sequence of bonded parti- 
cles is a self-assembled rod. At fixed temperature, the 
average rod length increases as the density increases and 
the polydisperse rods will undergo an nematic ordering 
transition [2?| . 

Since each site state is characterized by a three-state 
variable, we can rewrite Hamiltonian (TT]) in terms of new 
variables Si — 0, ±1, where where Si ± 1 represent the 
vertical (<?j = y) and the horizontal (er^ = x) orienta- 
tions, while Si = represents the empty state. Then, 
Hamiltonian |T]) reads 
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FIG. 2: Size dependence of the order parameter as a function 
of temperature for 9 = 0.525 and 8 = 1 (inset). 



This Hamiltonian has the same energy spectrum as (fT). 
Notice that the transformation Si — > —Si is not a sym- 
metry of the Hamiltonian ([2]) , since it is equivalent to a 
90° rotation of the lattice, but it is a symmetry of the sys- 
tem, since it left the partition function unchanged. The 
total number of adsorbed particles can be written 



(3) 



When N = M, we have Sf = 1 and the Hamiltonian 
3) results 



H({Si}) 



4 ^ 

<i,j> 



[1 + SiSj + (Si + Sj)(y.i*ij - x.r i:j )] 

(4) 

The last sum in Eq.(4) vanishes and therefore Hi = 
— T S<i j> SiSj+ constant. Hence, in that limit the 
present model reduces to the Ising one with coupling con- 
stant w Ismg = w/4. 



III. MONTE CARLO SIMULATIONS 



A. Monte Carlo Method 



scaling techniques [48[. The procedure is as follows. 
Starting with a random initial configuration (sites oc- 
cupied with concentration 9 = N/M and particle axis 
orientation chosen with probability 1/2), successive con- 
figurations are generated by attempting to move sin- 
gle particles (monomers). One of the two (translation 
or rotation) moves is chosen at random. In a trans- 
lation move, an occupied site and an empty site are 
randomly selected and their coordinates are established. 
Then, an attempt is made to interchange its occup ancy 
state with probability given by the Metropolis rule [49j : 
P = min {1, exp (— (3 AH)}, where AH is the difference 
between the Hamiltonians of the final and initial states 
and j3 = X/ksT (being ks the Boltzmann constant). For 
a rotation move, the rotational state of the chosen parti- 
cle (horizontal or vertical) is changed with a probability 
determined by Metropolis criteria. 

A Monte Carlo step (MCS) is achieved when 9 x M 
sites have been tested to change its occupancy state. 
Typically, the equilibrium state can be well reproduced 
after discarding the first 5 x 10 6 MCS. Then, the next 
6 x 10 8 MCS are used to compute averages. All calcula- 
tions were carried out using the parallel cluster BACO of 
Universidad Nacional de San Luis, Argentina. This fa- 
cility consists of 60 PCs each with a 3.0 GHz Pcntium-4 
processor and 90 PCs each with a 2.4 GHz Intel Core 2 
Quad processor. 

In order to follow the formation of the nematic phase 
from the isotropic phase, we use the order parameter de- 
fined in Ref. [H, 



\N V -N h \ 

Ny + N h 



(5) 



where Nf l (N v ) is the number of monomers aligned along 
the horizontal (vertical) direction (N = Nh + N v ). 

In our Monte Carlo simulations, we set the density 
6, varied the temperature T and monitored the order 
parameter S, which can be calculated as simple averages. 
The reduced fourth-order cumulant, Ul, introduced by 
Binder [43, was calculated as: 



U L (T) = 1 



(S 4 



3(5 2 ) ; 



(6) 



We have used a standard importance sampling MC 
method in the canonical ensemble |47l and finite-size 



where the thermal average (...), in all the quantities, 
means the time average throughout the MC simulation. 
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FIG. 3: Curves of U L vs T* for 6 = 0.525 (a) and 9 = 1 (b). 
From their intersections one obtained T* . In the insets, the 
data are plotted over a wider range of temperatures. 



B. Computational results 

The critical behavior of the present model has been 
investigated by means of the computational scheme de- 
scribed in the previous section and finite-size scaling anal- 
ysis mut. 

We start with the calculation of the order parameter 
plotted versus the reduced temperature T* = ksT/w for 
several lattice sizes (L — 60, 80, 100 and 120) and two 
values of coverage [9 = 0.525 [50|, Fig. 2 and 9 = 1, 
inset of Fig. 2]. As it can be observed, 8 appears as a 
proper order parameter to elucidate the phase transition. 
When the system is disordered (T* > T*, being T* the 
critical temperature) , all orientations are equivalents and 
8 is zero. In the critical regime (T* < T*), the particles 
align along one direction and 8 is different from zero. 

Hereafter we discuss the behavior of the critical tem- 
perature as a function of coverage. The standard theory 
of finite-size scaling allows for various efficient routes to 
estimate T* from MC data [47|, El]. One of these meth- 
ods, which will be used in this case, is from the tempera- 
ture dependence of Ul(T*), which is independent of the 
system size for T* = T* . In other words, T* can be found 
from the intersection of the curve Ul(T*) for different 
values of L, since U* = Ul{T*) =const. As an example, 
Fig. 3 shows the reduced four-order cumulants Ul plot- 
ted versus T* for the cases studied in Fig. 2. The values 
obtained for the critical temperature were T* — 0.250(2) 



(corresponding to 6 = 0.525) and T* = 0.567(2) (corre- 
sponding to 9 = 1). The procedure was repeated for 9 
ranging between and 1. The results, which are collected 
in Fig. 4 (a), represent the temperature-coverage phase 
diagram of the system. The critical line (squares and line 
in the figure) separates regions of isotropic and nematic 
stability. The different phases are shown schematically 
in parts (b), (c) and (d) of Fig. 4. 

With respect to the numerical results obtained by 
Tavares et al. at 9 = 0.2 and 9 — 0.4 [denoted with solid 
circles in Fig. 4 (a)], the agreement with the present data 
is very good. 

As it is well-known, the behavior of the reduced fourth- 
order cumulant as a function of temperature not only 
provides an accurate estimation of the critical tempera- 
ture T c in the infinite system, but also allows to make 
a preliminary identification of the order and univers ality 
class of the phase transition occurring in the system (47| . 
In the case of Fig. 3, and as it is shown in the insets, the 
curves exhibit the typical behavior of the cumulants in 
the presence of a continuous phase transition. Namely, 
the order parameter cumulant shows a smooth drop from 
2/3 to 0, instead of a characteristic deep (negative) min- 
imum, as in a first-order phase transition [47J. 

With respect to the value of the intersection point U* , 
two different behaviors can be visualized from Fig. 3. On 
one hand, at 9 = 0.525, the value obtained for U* (U* = 
0.639(5)) is consistent with the q = 1 Potts universality 
class [51| observed in Ref. [29|, where the system was 
studied at a fixed temperature (T* = 0.25). On the other 
hand, and as it is expected for 9 = 1, the fixed value 
of the cumulants, U* = 0.611(1), is consistent with the 
extremely precise transfer matrix calculation of U* = 
0.6106901(5) for the 2D Ising model Even though 

the value of U* may be taken as a first indication of 
universality, a detailed calculation of critical exponents is 
required for an accurate determination of the universality 
class along the critical line in Fig. 4 (a), and this will be 
subject of future research. 

Finally, Fig. 5 shows the nematic order parameter 8 
as a function of the coverage. The data correspond to 
T* = 0.25 and L = 100 As the density is increased 
above a critical value, the particles align along one di- 
rection and 8 increases continuously to one, remaining 
constant up to full coverage. In other words, nematic 
order survives until 9 = 1. This finding (1) allows us to 
discard the existence of a reentrant nematic transition at 
high densities as speculated in Ref. [13] and (2) indicates 
a substantial difference between the present system and 
that of monodisperse rigid rods without self-assembly, 
where a second nematic to isotropic phase transition is 
observed at high densities [H, [HJ • 
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FIG. 4: (a) Phase diagram of the model: our simulation data (squares and line) and additional points (circles) obtained from 
Monte Carlo simulation, carried out by Tavares et al. p2?| . (b) Schematic representation of the low-density nematic phase 
(point A in the figure), (c) Same as (b) for the intermediate-density disordered phase (point B in the figure), (d) Same as (b) 
for the high-density nematic phase (point C in the figure). 



IV. ANALYTICAL APPROXIMATIONS AND 
COMPARISON BETWEEN SIMULATED AND 
THEORETICAL RESULTS 

In this section we calculate the phase diagram within 
mean field and real space renormalization group ap- 
proaches. Let 



9 = ^I>-> 



respectively, where 
ensemble average. 



(9) 



means here a grand canonical 



/ 



1 

M/3 



■In 



(J) 



the grand canonical free energy, where H' = H — /j,N, 
H and N are given by Eqs. ^ and ©, and \i is the 
chemical potential. The orientational order parameter 
and coverage are then given by 



A. Mean-Field Approximation 

To obtain a mean field free energy $ for this problem 
we use the variational method [5|| , based on Bogoliubov 
inequality 



1 

M 



(8) 



f<<S> = fo + jj(H'-H' ) . 



(10) 



and 



where Hq is a trial Hamiltonian containing variational 
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FIG. 5: Nematic order parameter if as a function of the cov- 
erage. The data correspond to T* = 0.25 and L — 100. 
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FIG. 6: T* vs. (j,/w mean field phase diagram. TCP is a 
tricritical point. 
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We choose 



where rj is an effective field that breaks the orientational 
symmetry. Then 



> 2 ) - - In {1 + 2 coshQ3r?)} , 

P 

(11) 



where 



_ 2e^sinh(^) 
Wo i + 2 e^cosh(^)' [ ' 



and 



, 2 2e^cosh(/3 7? ) 

\^/o i + 2e ^ cosh(/377)' 1 j 

Minimizing Eq. ([TT|) we obtain the self consistent equa- 
tion 



V 



1 



-S 2 



(14) 



We see that the isotropic state rj = (S = 0) is a solution 
of Eq. (fl"4|) . At low temperatures Eq. (|14[) also presents 
ordered (nematic) solutions 77 7^ 0. Making a Landau 
expansion of Eq. (fT4"]) we obtain the following results: 



• There is a tricritical point at T t * =3/4 and /x t 

a: 

coverage at this point is dt = 1/2 



^pln2, where a 2 = 04 = and a 6 > 0. The 



• When /J, > nt there is a second-order transition line 
(a 2 = 0, 04 > 0) at 




(15) 



Along the critical line we have T* = 9(2 — 8). When 
/1 — > 00 we have — >• 1 and T* — > 1. Moreover, 
from Eqs. ([T2"T) and (fT4|) we obtain in this limit 
(5 = tanh(/3u>(5), i.e., the mean field equation for an 
Ising model, as expected. 

• When fj, < —w only the isotropic solution remains. 
It is easy to see that there is a level crossing at 
T* = between the empty state 6 = rj = and the 
completely ordered one 5 = X] = 1. 

• When —w<fi<fi t we have a first-order transition 
line (04 < and a 6 > 0) which can be calculated 
numerically by a Maxwell construction. 

In Fig. 6 we show the mean field phase diagram at the 
(/x, T*) space, which is qualitatively similar to that of the 
isotropic Blume- Emery- Griffiths (BEG) model [56|. The 
corresponding phase diagram in (9,T*) space presents 
a coexistence region between a low-coverage isotropic 
phase and a high-coverage nematic one at low tempera- 
tures. The presence of this coexistence region (first-order 
phase transition) is completely at variance with the ob- 
served numerical simulation results. 



B. Real Space Renormalization Group approach 

In order to obtain a more accurate analytical predic- 
tion for the phase diagram we apply the Real Space 
Renormalization Group (RSRG) scheme introduced by 
Niemeijer and van Leeuwen |1J, using four spin Kadanoff 
blocks and a double majority rule RG projection matrix. 
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The details of the RSRG implementation are given in 
the supplementary material |58[. The application of a 
truncation scheme allowed us to restrict the proliferation 
of interactions. Under this framework, closed recursion 
RG relations can be obtained for the the more general 
Hamiltonian compatible with the basic symmetry of the 
system, namely, a 90 degrees rotation of the lattice when 
Si —Si, that is 



isotropic 



n RG = hJ2sf+ [LSiSi+M SfS]] 

i <i,j> 

+ J2i u ( s i s j + s * s i) (y- f ij - ^)]( 16 ) 

<i,j> 

where H = ~(3H and h = /3fi. For U = the Hamil- 
tonian (j2"4")l corresponds to the BEG model [56]. For 
L = M = U = fiw/A we recover the model ([2]). 

The RG flow starting from the subspace (L, M, U, h) = 
(K/4:,K/4:,K/4,h), with K = (3w, is governed by the 
following fixed points. (i) Two attractors at I± = 
(0,0,0,±oo). They represent the high {(Sf) « 1) 
and the low ((Sf) <C 1) density isotropic phases re- 
spectively, (ti) One semi- unstable fixed point T\ = 
(0, 0, 0, — In 2). It is the locus of a surface in the 
(L,M,U,h) space that corresponds to a smooth contin- 
uation at high temperatures between both phases. (Hi) 
A line of attractive fixed points at (+oo, 0, 0, +oo). It 
is the locus of the ferromagnetic phase in the whole 
(L,M,U,h) space and we call it the N attractor. (iv) 
One non-trivial fixed point C% — (L c , 0, 0, +oo) with 

2V2- 



L 



= \ In 



1 



0.518612. It is the 



VlO + 572 

locus of a critical surface and corresponds to the critical 
point of the Ising model in the square lattice under the 
present approximation. The associated critical exponent 
results v = 1.0013..., in excellent agreement with the ex- 
act result v = 1 . The details of this analysis are given in 
the supplementary material [HI]. 

The phase diagram in the (K, h) space, obtained 
from the RG flow starting with (L,M,U,h) — 
(K/4, K/4, K/4, h), is shown in Fig. 7. We found a single 
critical line separating the nematic and isotropic phases 
(black continuous line in Fig. 7), which is in the basin 
of attraction of the fixed point C\. The nematic phase 
is in the basin of of attraction of N, while the isotropic 
phase is attracted either by I + or by i_. Points along 
the grey dashed line in Fig. 7 are attracted by the trivial 
fixed point T\, thus corresponding to a smooth continu- 
ation from low to high density isotropic phases, without 
phase transition. This line converges asymptotically to 
the critical line when h — > —00. Therefore, according to 
the present RG prediction the transition is second order 
for any finite temperature and it is in the universality 
class of the Ising model. The corresponding phase dia- 
gram in the {n/w,T*) = (h/K,l/K) space is shown in 
the inset of Fig. 7. 

Finally, to calculate the phase diagram in the (9,T*) 
space we computed numerically the coverage 9(h) along 
the critical line K c = K c (h) of Fig. 7. The results are 




nematic 



isotropic 



-10 



-5 
h 



FIG. 7: RG phase diagram in the (K, h) space. The black 
continuous line is attracted by the fixed point C\ and there- 
fore corresponds to a secondorder critical one. The grey 
dashed line is attracted by the fixed point T\ and corresponds 
to a smooth continuation between the high and low density 
isotropic phases. The inset shows the corresponding phase 
diagram in the (fi/w,T*) — (h/K, 1/K) space. 



presented in Fig. 8 and the details of the calculation are 
given in the supplementary material [58|. 



Comparison Between Theoretical and 
Simulated Results 



In Fig. 8 we compare the critical lines obtained by MC 
(open squares joined by lines) and RSRG (solid circles), 
together with the analytical approximation developed by 
Tavares et al. [27| (solid line). While qualitatively sim- 
ilar to the MC result, we see that the present RSRG 
approximation systematically underestimates the critical 
temperature. Concerning the comparison with Tavares et 
al. results 27], quantitative and qualitative differences 
have been found between the analytical and the simula- 
tion data. In fact, the theory overestimates the critical 
temperature in all range of coverage, confirming the pre- 
dictions in Ref. [27J . For small values of 9, small differ- 
ences appear between simulation and theoretical results; 
however, the disagreement turns out to be significantly 
large for larger 9 , s. 

In the particular case of 9 = 1, the Tavares et al. the- 
ory predicts a critical temperature of T* = [ln(3/2)] _1 w 
2.466, whereas the value calculated by MC simulations 
is T* — 0.567(2). These results can be compared with 
the exact value of the critical temperature at full cover- 
age T* = - [2\n(V2- l)] _1 fa 0.567, (see section EJ). 
This result is consistent with that calculated by MC sim- 
ulations, which reinforces the robustness of the present 
computational scheme. 
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e 

FIG. 8: Comparison between numerical and theoretical esti- 
mates of the phase diagram in the (9, T*) phase diagram. 

V. CONCLUSIONS 

In summary, we have addressed the temperature- 
coverage phase diagram of self-assembled rigid rods on 
square lattices. By using Monte Carlo simulations, mean- 
field theory and a renormalization group approach, we 
obtained and characterized the critical line which sepa- 
rates regions of isotropic and nematic stability. Several 
conclusions can be drawn from the present results. 

First, a simulation test of the theory developed by 
Tavares et al. [27| was carried out. The results showed 
that the theory overestimates the critical temperature 
in all range of coverage, confirming the predictions in 
Ref. [27]. For small values of 9, small differences ap- 
pear between simulation and theoretical results; however, 
the disagreement turns out to be significantly large for 
larger #'s. On the other hand, the RSRG approach re- 
produces qualitatively the shape of the critical line, but 
systematically underestimates the critical temperature. 
Concerning this last calculation, the main prediction is 
that the critical properties of the whole line are associ- 
ated to a unique second-order fixed point, confirming the 
continuous nature of the transition. However, it must 



be pointed out that it predicts that the whole line is in 
the universality class of the d — 2 ferromagnetic Ising 
model, at variance with Monte Carlo numerical calcula- 
tions predicting that the transition at 9 « 1/2 belongs to 
the q — 1 Potts universality class [29( . While the present 
RSRG results are not conclusive, due to the approximate 
character of the approach, they indicate that further re- 
search is required to clarify this point. 

On the other hand, the behavior of the order param- 
eter allowed to discard the existence of a reentrant ne- 
matic transition at high densities as speculated in Ref. 
[27| . This finding indicates a substantial difference be- 
tween the present system and that of monodisperse rigid 
rods without self-assembly, where a second nematic to 
isotropic phase transition is observed at high densities 
[Mill 

Concerning the MF results, the prediction of a first- 
order transition line and a tricritical point is not sur- 
prising, due to the close relationship between the present 
model and the BEG one, as evidenced by the Eq. ([2|). In- 
deed, the generalized form contains both first-order 
and tricritical fixed points, but the RSRG results show 
that in d = 2 the anisotropic character of the interac- 
tions drive the RG flow of the present system outside 
their basins of attraction. However, in three dimensional 
systems the IN transition is usually first-order [401 ] . On 
the other hand, from the exact mapping into the isotropic 
Ising model at full coverage one could expect a second- 
order transition for high values of the coverage, even in 
three dimensions. Hence, the MF prediction of a tricrit- 
ical point is probably correct for d > 2. 
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RSRG Method 

Consider the Kadanoff blocks of size Nt = b 2 = 4 
shown in Fig. 9. Let's denote by S' T the block spin as- 
sociated to the block I and sj the set of lattice spins 
belonging to the block /: S] = {si} with i € I. Let's 
also denote by S' and s the complete sets of block and 
lattice spins respectively. We can express H = Ho + V, 
where Ho = J2i^-i( s i) contains all the interactions be- 
tween spins belonging to the block I and V all the in- 
terblock interactions. Introducing an RG projection ma- 
trix P(S', s) — YIt Pi(S'j, Si), an average of an arbitrary 
function X(S',s) as 
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FIG. 9: Kadanoff blocks of size b — 2 for the square lattice. 



(X) (S>) = -L £ P(S', S )e H ^X(S>, s) (17) 
^o — 

s 

where Zo = Yii %o with 



The simplest RG approach within Niemeijer and van 
Leeuwen [l| scheme consist into the identification 



L' 
M' 
U' 
hi 



2La( 

2Ma\ 

2U a\a2 

8M a,2 as + 04, 



(25) 
(26) 
(27) 
(28) 



together with 



H'(S')+C = \nZ + (V) . 



(18) 



where H'(S') is the block Hamiltonian and C is a spin- 
independent constant. This uncontrolled approximation 
results from the truncation to the first-order cumulant [l| 
of (exp(V)) . Using then of the double majority rule 
RG projection matrix Pj(S'j,si) introduced by Berker 
and Wortis for the pure isotropic Blume-Emery-Griffiths 
(BEG) model [2], it is easy to see that (Su) = aiS'j, 
(^ii)q = a 2 S? + a 3> an d m -^o = a iSf + »5, where 



ai = (S , i j> | Sf , =1 

«2 = (Su) Q \ s , i=1 - (S?/> | s , =0 

°3 = ( S ?i)o\s>=o 

a 4 = lnZ£| s}=1 - ln^| s , =0 

ls;=o • 



a 5 = In Z\ I 



(19) 
(20) 
(21) 
(22) 
(23) 



Applying this scheme to the Hamiltonian 



H RG = hJ2sf+ J2 [LSiSi + MS?S] 

i <ij> 

+ U (S^S, + S 2 S. t ) (y.r l3 - x.ni)] , (24) 



Defining 



we obtain 



a-3 
and 



g = C/N=(Ma 2 3 +a 5 /A). 



Bi(L, M, U, h) 
B 2 (L,M,U,h) 



Zj v\s' I =o 

Z o\s'=l 



(29) 



26^ + 26 



2h 



2e 



2H+M-L 



+ 2e 



2h+M+L 



cosh(2[/) 



Bi(L,M,f7,ft) 



B 2 (L,M,U,h) 
n B! (L,M,U,h) 
In 5i (L, M,C/,/i) 
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B 2 (L,M,U,h) 



- e 2h + 2 e 4(/l+M) + e 3'i+2(M-L) , g £ 3h+2(M+L) , g 4(h+M+L)_ 



+ e 



2/i+M+L 



cosh(2Z7) + 2 e 



3h+2M 



cosh(2C/)] 



« 2 



03 = 



B 2 (L,M,U,h) 



e 2h + g e 4(M+^) + e 2/ l+ M-L + g e 3/.+2M cosh(2i) + 2 ^(M+h) cosh(4L) + 



e 



2h+M+L 



cosh(2C7) + 6e 



2M+3A 



sh(2£/)] -a 3 (L,M,U,h) 



2e ll + 2e 



2 1 1 



2e 



2fi+M-L 



2e 



2/i+M+L 



cosh(2L/) 



Bi 
B 2 



1 + 8 e h + 4 e 



„2/i 



4e 



2/i+M-L 



4e 



2h+M+L 



cosh(2£7) 



2 e 2ft + 6 e 4 ( M+,i ) + 2 e 2h + M - £ + 8 e 3 ' i+2M cosh(2L) + 2 e 4 ( M +^ cosh(4L) 
2e 2h+M+L cosh(2[/) + 8e 3/ l+ 2Af cosh{ 2U). 

I 



RSRG flow and flxed points structure 

We found that all the relevant fixed points of the re- 
cursion equations lie in the BEG subspace U = 0. The 
RG flow and the fixed points structure in the U = sub- 
space is qualitatively similar to that obtained in Ref . Q , 
including first and second-order surfaces, as well as tri- 
critical and critical endpoint lines • We focused only on 
those fixed points relevant to the present problem namely, 
those which govern the RG flow starting from the sub- 
space (L,M,U,h) = (K/4,K/4:,K/4,h), with K = /3w. 
The whole flow starting from that subspace is attracted 
by two subspaces invariant under RG: L = M = U = 
and (M, h) = (0,+co). 



Flow in the L = M = U — subspace 

The recursion relations in this case reduce to b! = 
04(0, 0, 0, h). This RG equation presents three fixed 
points: two attractors at h — ±00, which are the loci 
of the high (h — 00) and low (h = — 00) density isotropic 
phases respectively and one unstable high temperature 
fixed point at h = — In 2. The first two fixed points are 
attractors in the complete (L, M, U, h) space and we will 
call them 7+ and 7_. They represent the high ((Sf) ~ 1) 
and the low ((Sf) 1) density isotropic phases respec- 
tively. The fixed point T\ = (0, 0, 0, — In 2) is the locus of 
a surface in the [L, M, U, h) space that corresponds to a 
smooth continuation at high temperatures between both 
phases. 



Flow in the (M, h) = (0, +00) subspace 

This subspace corresponds to an anisotropic Ising 
model, since in this limit the Si = state has zero prob- 
ability. The recursion relations reduce in this case to 



V = 2Ld(L) 2 
U' = Ud{L) 



with 



„4L 



d(L) = lim ai(L,0,U,h) 

h— >oo 



6 + 2 cosh(4L) ' 



(30) 
(31) 



(32) 



Since L' = L'(L), independently of the parameter U, the 
whole flow is governed by the RG equation correspond- 
ing to the isotropic Ising model. This equation has a 
non trivial fixed point at d(L) = 1/V2, whose solution 

VTo~ 



is L r 



1 



0.518612, corre- 



2\/2+ V10 + 5V2 

sponding to the critical point of the Ising model in the 
square lattice under the present approximation (compare 
with exact Onsager result L c = \ ln(l + w 0.44069). 
We will call this fixed point C\. The critical exponent v 
is given by 



In b 



A = 



We obtain v = 1.0013..., in excellent agreement with the 
exact result v = 1. The RG recursion equation also 
has two attractors: J+ (L — 0) and the isotropic fer- 
romagnetic fixed point L — 00 (T — 0). At L = L c we 
have another invariant line at the (L, U) space, whose 
RG equation is U' = U' j\p2. This recursion relation has 
only trivial fixed points: one attractor at U = and one 
unstable at U — +00. The line L = is also invariant 
and have the same fixed points. Finally, we have that 
liirt£_>. 00 d(L) = 1 Hence, U' = U and the whole line 
L = +00 is a line of fixed points. This is the locus of the 
ferromagnetic phase in the whole (L, M, U, h) space and 
we will call it the N attractor. In Fig. 10 we show the 
flow diagram in the complete (U, L) space. 
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u 

FIG. 10: RG flow in the (L, 0, U, +oo) invariant subspace. 

RSRG Coverage calculation 

The coverage can be expressed as 



9(K, h) = -13 



df(K,h) 
dh ' 



(33) 



Let K = (L, M, U, h) be the parameters vector of Hamil- 



tonian (|24[) . From the renormalization group transforma- 
tion we have the followingrelation after n applications 
of the RG transformation [1| 



f(K 



i " 



md 



g(K n 



-nd 



f(K n ) 



(34) 



where K m is the parameters vector after m applications 
of the RG transformation, Kq is the initial value and 
g(K) = C/N is given by Eq. ([2"5]) . Since 9 is not singular 
at the critical line, we can assume that the singular part 
of the free energy will make no contribution to Eq. (|3"4")l 
and therefore the derivative of the second term in the 
right hand of the previous expression vanishes when n — > 
oo. Therefore, we can express 



0(K, h) 



dh 



.m=0 



%d g(K m ) 



K = (K/4,K/4,K,4,h) 



(35) 

Computing numerically the above sum and taking the 
numerical derivative we obtain the critical line T* vs. 9 
shown in Fig. 8 of the manuscript. 
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